#include "navierstokes.hpp"

// variational formulation language
//#include <life/lifevf/vf.hpp>

#include <ns_cl.cpp>

namespace Feel
{
namespace FSI
{
    template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
    void
    NavierStokes<Dim, Order, Entity>::updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J, vector_ptrtype& R)
    {

        using namespace Feel::vf;

        if (M_LookComment) std::cout <<"[NavierStokes] : updateJacobian start \n";
        boost::timer btime;
        btime.restart();

        this->changeRepository();

        //mesh_ptrtype mesh = Xh->mesh();
        auto mesh = this->mesh();
        auto Xh = this->functionSpace();


        element_fluid_type U( Xh, "u" ); U = *X;
        element_fluid_type V( Xh, "v" );
        element_fluid_0_type u = U.template element<0>();
        element_fluid_0_type v = V.template element<0>();
        element_fluid_1_type p = U.template element<1>();
        element_fluid_1_type q = V.template element<1>();

        //--------------------------------------------------------------------------------------------------//
#if 1
        //Structure data
        double EE = 21*1e5;
        double sigma = 0.28;
        // !!!! Connerie avec le mu !!!!!!!!!!!!!!!!!! 
        double mu = EE/(2*(1+sigma));
        double lambda = EE*sigma/((1+sigma)*(1-2*sigma));

        auto Iddd = oneX()*trans(oneX()) + oneY()*trans(oneY());
        auto Fv = Iddd + gradv(*M_aux);
        auto Ev = 0.5*(gradv(*M_aux)+trans(gradv(*M_aux)) ) + 0.5*trans(gradv(*M_aux))*gradv(*M_aux);
        auto Sv = lambda*trace(Ev)*Iddd + 2*mu*Ev;
#endif
        //--------------------------------------------------------------------------------------------------//

        //Tenseur des deformations
        auto deft = 0.5*(gradt(u)+trans(gradt(u)) );
        auto def = 0.5*(grad(v)+trans(grad(v)) );
        //Matrice Identitee
        auto Id = oneX()*trans(oneX()) + oneY()*trans(oneY());
        // Tenseur des contraintes (trial)
        auto Sigmat = -idt(p)*Id + 2*M_mu*deft;
        // Tenseur des contraintes (test)
        auto Sigma = -id(q)*Id + 2*M_mu*def;
        //force volumique :
        //auto f = -rho*cst(9.81)*vec(cst(0.),cst(1.));
        // Vitesse d'entree du fluide
        //auto Vit_Entree = 10*Vitesse_Entree*( Py()*(1.0-Py()) )*oneX();

        //std::vector<void*> listTemp(1);
        //auto ola = (void *)(& Vit_Entree);
        //listTemp[0]=ola;

        //--------------------------------------------------------------------------------------------------//
        //--------------------------------------------------------------------------------------------------//
        //--------------------------------------------------------------------------------------------------//

        form2( Xh,Xh, J , _init=true );

        //terme non lineaire
        bool M_isNLTerme=true;
        if (M_isNLTerme)
            {
                form2( Xh,Xh, J ) +=
                    integrate (elements(mesh), _Q<3*Order>(),
                               rho*trans(id(v))*(gradv(u))*idt(u));
                form2( Xh,Xh, J )  +=
                    integrate (elements(mesh), _Q<3*Order>(),
                               rho*trans(id(v))*(gradt(u)*idv(u)) );

                if (M_isMoveDomain)
                    {
                        form2( Xh,Xh, J )  +=
                            integrate (elements(mesh), _Q<3*Order>(),
                                       -rho*trans(idv(*M_mesh_velocity))*(gradt(u))*id(v) );
                    }
            }

#if 1
        // sigma : grad(v) sur Omega
        form2( Xh, Xh, J ) +=
            integrate( elements(mesh), _Q<2*(Order)>(),
                       + trace(Sigmat*trans(grad(v))) );
#if 0
        form2( Xh, Xh, J ) +=
            integrate( elements(mesh),// _Q<2*(Order)>(),
                       - rho*divv(*M_mesh_velocity)*trans(idt(u))*id(v) );
#endif
#else
        form2( Xh, Xh, J ) +=
            integrate( elements(mesh),
                       + 2*M_mu*trace(deft*trans(grad(v))) );
        form2( Xh, Xh, J ) +=
            integrate( elements(mesh),
                       + div(v)*idt(p) );


#endif


        // terme assurant l'incompressibilité
        form2( Xh, Xh, J ) +=
            integrate( elements(mesh), _Q<2*(Order-1)>(),
                       + divt(u)*id(q) );

        //--------------------------------------------------------------------------------------------------//

        auto clDef = CL_DEF_NS(M_time);

        if (M_weakCL)
            {
                ForEachCL( clDef, cl::dirichlet_vec,
                           form2( Xh, Xh, J ) +=
                           integrate( markedfaces(mesh,mesh->markerName(PhysicalName)), //_Q<2*Order>(),
                                      - trans(Sigmat*N())*id(v)
                                      + penalbc*trans(idt(u))*id(v)/hFace() ) );

                ForEachCL( clDef, cl::paroi_mobile,
                           form2( Xh, Xh, J ) +=
                           integrate( markedfaces(mesh,mesh->markerName(PhysicalName)), //_Q<2*Order>(),
                                      - trans(Sigmat*N())*id(v)
                                      + penalbc*trans(idt(u))*id(v)/hFace() ) );
#if 0
                // terme de bord classique
                form2( Xh, Xh, J ) +=
                    integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                               - trans(Sigmat*N())*id(v) );
                if(M_isFixeBoundary)
                    form2( Xh, Xh, J ) +=
                        integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
                                   - trans(Sigmat*N())*id(v) );

                if (!M_isMoveDomain)
                    {
                        form2( Xh, Xh, J ) +=
                            integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                                       - trans(Sigmat*N())*id(v) );
                        form2( Xh, Xh, J ) +=
                            integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                                       - trans(Sigmat*N())*id(v) );
                    }
                //terme de bord weak boundary
                form2( Xh, Xh, J ) +=
                    integrate( markedfaces(mesh,mesh->markerName("Entree")), _Q<2*Order>(),
                               - trans(Sigma*N())*idt(u)
                               + penalbc*trans(idt(u))*id(v)/hFace() );
                if(M_isFixeBoundary)
                    form2( Xh, Xh, J ) +=
                        integrate( markedfaces(mesh,mesh->markerName("Paroi")), _Q<2*Order>(),
                                   - trans(Sigma*N())*idt(u)
                                   + penalbc*trans(idt(u))*id(v)/hFace() );
                form2( Xh, Xh, J ) +=
                    integrate( markedfaces(mesh,mesh->markerName("ParoiH")), _Q<2*Order>(),
                               //- trans(Sigma*N())*idt(u)
                               + penalbc*trans(idt(u))*id(v)/hFace() );
                form2( Xh, Xh, J ) +=
                    integrate( markedfaces(mesh,mesh->markerName("ParoiB")), _Q<2*Order>(),
                               //- trans(Sigma*N())*idt(u)
                               + penalbc*trans(idt(u))*id(v)/hFace() );
#endif

            }
        // terme necessaire pour condition forte avec petsc
        form2( Xh, Xh, J ) +=
            integrate( elements( mesh ),// _Q<3*uOrder-1+3*(GeoOrder-1)>(),
                       +0*div(v)*idt(p)+0*divt(u)*id(q) +0*idt(p)*id(q) );

        //bool M_isStat=true;
        if (!M_isStat)
            {
                form2( Xh, Xh, J ) +=
                    integrate( elements(mesh), _Q< 2*Order >(),
                               + rho*trans(idt(u))*id(v)*M_bdf_fluid->polyDerivCoefficient(0) );
            }

#if 0
        int stabilisation = this->application()->vm()["stab"].template as<int>();
        if (stabilisation==1) {
            //terme de stabilisation
            form2( Xh, Xh, J ) +=
                integrate( internalfaces(mesh), _Q< 3*Order >(),
                           + 0.025*vf::pow(hFace(),2.0)*( trans(idt(u))*N() )*(trans(jump(grad(u)))*jump(grad(v)) )  );
            //+ 0.025*vf::pow(hFace(),2.0)*( abs( trans(idv(u))*N() ))*(trans(jumpt(gradt(u)))*jump(grad(v)) )  );
            //            + (0.05+vf::pow(hFace(),2.0))*( (trans(idv(u))*N())*( trans(idv(u))*N() ))*(trans(jumpt(gradt(u)))*jump(grad(v)) )  );
            //( (trans(idt(u))*N())*( trans(id(u))*N() )) );

            form2( Xh, Xh, J ) +=
                integrate( elements(mesh), _Q< 3*Order >(),
                           + 0.000*vf::pow(hFace(),2.0)*( abs( trans(idv(u))*N() ))*(trans(jumpt(gradt(u)))*jump(grad(v)) )  );
        }
#endif


        J->close();

        if (!M_weakCL)
            {
                auto dirElem = this->getDirichletElement();
                ForEachCL( clDef,cl::paroi_mobile,
                           form2( Xh, Xh, J ) +=
                           on( markedfaces(Xh->mesh(), Xh->mesh()->markerName(PhysicalName)) , *dirElem, R, vf::vec(vf::cst(0.),vf::cst(0.)) ) );

                ForEachCL( clDef,cl::dirichlet_vec,
                           form2( Xh, Xh, J ) +=
                           on( markedfaces(Xh->mesh(), Xh->mesh()->markerName(PhysicalName)) , *dirElem, R, vf::vec(vf::cst(0.),vf::cst(0.)) ) );

            }

        double timeElaps = btime.elapsed();
        if (M_LookComment) std::cout <<"[NavierStokes] : updateJacobian finish in " << timeElaps << "s\n";

    }

} // end namespace NavierStokes

} // Feel


// instantiation
template class Feel::FSI::NavierStokes<2,2,Feel::Simplex>;
//template class Feel::NavierStokes<2,4,Feel::Simplex>;

